// Implementation of the age-structured model of pertussis transmission
// Model includes maternal immunization, primary vaccination (2nd age group, coded as 1), infant booster (3rd age group, coded as 2) and the option of a pre-school booster (age group coded as 7)

//start_globs
// Expit transform for parameters constrained in interval [a,b]
static double expitCons(double x, double a, double b) {
double out = (a + b * exp(x)) / (1.0 + exp(x));
if(ISNAN(out)) out = (b + a * exp(-x)) / (1.0 + exp(-x)); // If x=+Inf, must return b
return out;
}

// Logit transform for parameters constrained in interval [a,b]
static double logitCons(double x, double a, double b) {
x = (x <= a) ? a : (x >= b ? b : x);
double out = log((x - a) / (b - x));
return out;
}

// Probability of transition for event of rate r during time step delta_t
// p = 1 - exp(-r * delta_t)
static double pTrans(double r, double delta_t) {

// r: event (instantaneous rate)
// delta_t: time step
// Returns: probability of transition
double out = 1.0 - exp(-r * delta_t);
return out;
}
//end_globs

//start_dmeas
// Observation model: evaluation of log-likelihood
double logL = 0.0;
int i; 
double *CI1_vec = (double *) &CI1_1; // Cumulative no
double *CI2_vec = (double *) &CI2_1; // Cumulative no
double *CI3_vec = (double *) &CI3_1; // Cumulative no
double *CI4_vec = (double *) &CI4_1; // Cumulative no
double *CI5_vec = (double *) &CI5_1; // Cumulative no
double *CIobs_vec = (double *) &CIobs_1; // Observed no of cases
double *rho_vec = (double *) &rho_1; // Reporting probabilities
for(i = 0; i < nages_mod; i++) {
	logL += dnbinom_mu(nearbyint(CIobs_vec[i]), 1.0 / kC, rho_vec[i] * (CI1_vec[i] + CI2_vec[i] + CI3_vec[i] + CI4_vec[i] + CI5_vec[i]), 1);
}
lik = (give_log) ? (logL) : exp(logL);
//end_dmeas

//start_rmeas
// Observation model: generate observations
int i; 
double *CI1_vec = (double *) &CI1_1; // Cumulative no
double *CI2_vec = (double *) &CI2_1; // Cumulative no
double *CI3_vec = (double *) &CI3_1; // Cumulative no
double *CI4_vec = (double *) &CI4_1; // Cumulative no
double *CI5_vec = (double *) &CI5_1; // Cumulative no
double *CIobs_vec = (double *) &CIobs_1; // Observed no of cases
double *rho_vec = (double *) &rho_1; // Reporting probabilities
for(i = 0; i < nages_mod; i++) {
	CIobs_vec[i] = rnbinom_mu(1.0 / kC, rho_vec[i] * (CI1_vec[i] + CI2_vec[i] + CI3_vec[i] + CI4_vec[i] + CI5_vec[i]));
}
//end_rmeas

//start_skel
// Skeleton of deterministic model
int i, j, i2;

//Pointers to state variables
double *VM_vec = (double *) &VM_1;
double *VMbar_vec = (double *) &VMbar_1;
double *V_vec = (double *) &V_1;
double *M_vec = (double *) &M_1;
double *MV_vec = (double *) &MV_1;
double *SMbar_vec = (double *) &SMbar_1;
double *S_vec = (double *) &S_1;
double *E_vec = (double *) &E_1;
double *I_vec = (double *) &I_1;
double *SV_vec = (double *) &SV_1;
double *SVMbar_vec = (double *) &SVMbar_1;
double *SVM_vec = (double *) &SVM_1;
double *R_vec = (double *) &R_1;
double *CI1_vec = (double *) &CI1_1;
double *CI2_vec = (double *) &CI2_1;
double *CI3_vec = (double *) &CI3_1;
double *CI4_vec = (double *) &CI4_1;
double *CI5_vec = (double *) &CI5_1;

//Pointers to derivatives
double *DVM_vec = (double *) &DVM_1;
double *DVMbar_vec = (double *) &DVMbar_1;
double *DV_vec = (double *) &DV_1;
double *DM_vec = (double *) &DM_1;
double *DMV_vec = (double *) &DMV_1;
double *DSMbar_vec = (double *) &DSMbar_1;
double *DS_vec = (double *) &DS_1;
double *DE_vec = (double *) &DE_1;
double *DI_vec = (double *) &DI_1;
double *DSV_vec = (double *) &DSV_1;
double *DSVMbar_vec = (double *) &DSVMbar_1;
double *DSVM_vec = (double *) &DSVM_1;
double *DR_vec = (double *) &DR_1;
double *DCI1_vec = (double *) &DCI1_1;
double *DCI2_vec = (double *) &DCI2_1;
double *DCI3_vec = (double *) &DCI3_1;
double *DCI4_vec = (double *) &DCI4_1;
double *DCI5_vec = (double *) &DCI5_1;

// Vector of parameters
double *q_vec = (double *) &q_1; // age-specific clinical fraction
double *N_vec = (double *) &N_1; // age-specific population size
double *delta_vec = (double *) &delta_1; // age-specific aging rates 
double lambda_vec[nages_mod];

// Calculate vaccine coverages
double p0_t = (t < tM) ? 0.0 : v0; // Maternal immunization
double p1_t = (t < tV) ? 0.0 : v1; // Primary series
double p2_t = (t < tV) ? 0.0 : v2; // Booster doses

// Re-create contact matrix
double CM_mat[nages_cmat][nages_cmat]; 
double *CM_vec = (double *) &CM_1; // Contact matrix stacked as vector (line by line), size nages_cmat*nages_cmat
for(i = 0; i < nages_cmat; i++) {
	for(j = 0; j < nages_cmat; j++) {
		CM_mat[i][j] = CM_vec[nages_cmat * i + j];
	}
}

// Compute seasonality factors in 5-10y and 10-20 y
CM_mat[1][1] *= seasC; // 5--10 y
CM_mat[2][2] *= seasT; // 10--15 y
CM_mat[3][3] *= seasT; // 15--20 y

// Create augmented contact matrix
double CM_mat_aug[nages_mod][nages_cmat]; // "Augmented" contact matrix with 2 extra rows for infant classes
double I_tilde[nages_cmat], N_tilde[nages_cmat]; // No infected and population size without infant age classes
I_tilde[0] = I_vec[0] + I_vec[1] + I_vec[2] + I_vec[3] + I_vec[4] + I_vec[5]; // Sum no infected in 0-4 mo, 4-12 mo, and 1-4 y
N_tilde[0] = N_vec[0] + N_vec[1] + N_vec[2] + N_vec[3] + N_vec[4] + N_vec[5]; // Sum population in 0-4 mo, 4-12 mo, and 1-4 y
for(i = 1; i < nages_cmat; i++) {
	I_tilde[i] = N_tilde[i] = 0.0;
	for(j = 1; j <= 5; j++) {
		I_tilde[i] += I_vec[5 * i + j];
		N_tilde[i] += N_vec[5 * i + j];
	}
}

// Add rows for the augmented contact matrix
for(i = 0; i < nages_mod; i++) {
	for(j = 0; j < nages_cmat; j++){
		if(i <= 5) {
			CM_mat_aug[i][j] = CM_mat[0][j];
		} else {
			i2 = (int) floor((i - 1) / 5); // Range: 1--14
			CM_mat_aug[i][j] = CM_mat[i2][j];
		}
	}
}

// Calculate force of infection in every age group
for(i = 0; i < nages_mod; i++) {
	lambda_vec[i] = 0.0; // Initialize values
	for(j = 0; j < nages_cmat; j++){
		lambda_vec[i] += CM_mat_aug[i][j] * (I_tilde[j] + iota) / N_tilde[j];
	}
	lambda_vec[i] *= 1 * q_vec[i]; //adjust R0 here with 0.5 or 1.5 for sensitivity analyses
}

// Blunting 
double epsilon_bar = (1.0 - b1) * epsilon; // Blunted vaccine effectiveness
double alphaV_bar = alphaV / (1.0 - b2); // Blunted duration of vaccine protection

if(debug) {
	Rprintf("t=%.2f, seasC = %.2f, seasT = %.2f\n\n", t, seasC, seasT);
}

// Derivatives

// Equations in newborns
DVM_vec[0] = 0.0;
DVMbar_vec[0] = 0.0;
DV_vec[0] = 0.0;
DM_vec[0] = p0_t * epsilonM * mu * N_tot - (tau + delta_vec[0]) * M_vec[0];
DMV_vec[0] = 0.0;
DSMbar_vec[0] = p0_t * (1 - epsilonM) * mu * N_tot + tau * M_vec[0] - (lambda_vec[0] + delta_vec[0]) * SMbar_vec[0];
DS_vec[0] = (1.0 - p0_t) * mu * N_tot  - (lambda_vec[0] + delta_vec[0]) * S_vec[0];
DE_vec[0] = lambda_vec[0] * (S_vec[0] + SV_vec[0] + SMbar_vec[0] + SVMbar_vec[0] + SVM_vec[0]) - (sigma + delta_vec[0]) * E_vec[0];
DI_vec[0] = sigma * E_vec[0] - (gamma + delta_vec[0]) * I_vec[0];
DSV_vec[0] = 0.0;
DSVMbar_vec[0] = 0.0;
DSVM_vec[0] = 0.0;
DR_vec[0] = gamma * I_vec[0] - delta_vec[0] * R_vec[0];
DCI1_vec[0] = lambda_vec[0] * SMbar_vec[0];
DCI2_vec[0] = lambda_vec[0] * S_vec[0];
DCI3_vec[0] = lambda_vec[0] * SV_vec[0];
DCI4_vec[0] = lambda_vec[0] * SVMbar_vec[0];
DCI5_vec[0] = lambda_vec[0] * SVM_vec[0];

// Equations in older age groups
for(i = 1; i < nages_mod; i++) {
	DVM_vec[i] = delta_vec[i - 1] * VM_vec[i - 1] - (alphaV_bar + delta_vec[i]) * VM_vec[i]; 
	DVMbar_vec[i] = delta_vec[i - 1] * VMbar_vec[i - 1] - (alphaV + delta_vec[i]) * VMbar_vec[i];
	DV_vec[i] = delta_vec[i - 1] * V_vec[i - 1] - (alphaV + delta_vec[i]) * V_vec[i]; 
	DM_vec[i] = delta_vec[i - 1] * M_vec[i - 1] - (tau + delta_vec[i]) * M_vec[i];
	DMV_vec[i] = delta_vec[i - 1] * MV_vec[i - 1] - (tau + delta_vec[i]) * MV_vec[i];
	DSMbar_vec[i] = delta_vec[i - 1] * SMbar_vec[i - 1] + tau * M_vec[i] - (lambda_vec[i] + delta_vec[i]) * SMbar_vec[i];
	DS_vec[i] = delta_vec[i - 1] * S_vec[i - 1] - (lambda_vec[i] + delta_vec[i]) * S_vec[i];
	DE_vec[i] = delta_vec[i - 1] * E_vec[i - 1] + lambda_vec[i] * (S_vec[i] + SMbar_vec[i] + SV_vec[i] + SVMbar_vec[i] + SVM_vec[i]) - (sigma + delta_vec[i]) * E_vec[i];
	DI_vec[i] = delta_vec[i - 1] * I_vec[i - 1] + sigma * E_vec[i] - (gamma + delta_vec[i]) * I_vec[i];
	DSV_vec[i] = delta_vec[i - 1] * SV_vec[i - 1] + alphaV * V_vec[i]  - (lambda_vec[i] + delta_vec[i]) * SV_vec[i];
	DSVMbar_vec[i] = delta_vec[i - 1] * SVMbar_vec[i - 1] + alphaV * VMbar_vec[i] - (lambda_vec[i] + delta_vec[i]) * SVMbar_vec[i];
	DSVM_vec[i] = delta_vec[i - 1] * SVM_vec[i - 1] + alphaV_bar * VM_vec[i] + tau * MV_vec[i]  - (lambda_vec[i] + delta_vec[i]) * SVM_vec[i];
	DR_vec[i] = delta_vec[i - 1] * R_vec[i - 1] + gamma * I_vec[i] - delta_vec[i] * R_vec[i];
	DCI1_vec[i] = lambda_vec[i] * SMbar_vec[i];
	DCI2_vec[i] = lambda_vec[i] * S_vec[i];
	DCI3_vec[i] = lambda_vec[i] * SV_vec[i];
	DCI4_vec[i] = lambda_vec[i] * SVMbar_vec[i];
	DCI5_vec[i] = lambda_vec[i] * SVM_vec[i];
}

// Vaccination

// Primary course (2nd age group, i = 1)
DVM_vec[1] += p1_t * epsilon_bar * delta_vec[0]* M_vec[0];
DVMbar_vec[1] += p1_t * epsilon * delta_vec[0]* SMbar_vec[0];
DV_vec[1] += p1_t * epsilon * delta_vec[0] * S_vec[0]; 
DM_vec[1] -= p1_t * delta_vec[0] * M_vec[0];
DMV_vec[1] += p1_t * (1.0 - epsilon_bar) * delta_vec[0] * M_vec[0]; 
DSMbar_vec[1] -= p1_t * delta_vec[0] * SMbar_vec[0]; 
DS_vec[1] -= p1_t * delta_vec[0] * S_vec[0];  
DSV_vec[1] += p1_t * (1.0 - epsilon) * delta_vec[0] * S_vec[0];
DSVMbar_vec[1] += p1_t * (1.0 - epsilon) * delta_vec[0] * SMbar_vec[0];

// Pediatric booster (3rd age group, i = 2)
DVM_vec[2] += p2_t * epsilon_bar * delta_vec[1] * MV_vec[1] + p2_t * epsilon_bar * delta_vec[1] * SVM_vec[1];
DVMbar_vec[2] += p2_t * epsilon * delta_vec[1] * SVMbar_vec[1];
DV_vec[2] += p2_t * epsilon * delta_vec[1] * SV_vec[1]; 
DMV_vec[2] -= p2_t * epsilon_bar * delta_vec[1] * MV_vec[1];   
DSV_vec[2] -= p2_t * epsilon * delta_vec[1] * SV_vec[1];
DSVMbar_vec[2] -= p2_t * epsilon * delta_vec[1] * SVMbar_vec[1];
DSVM_vec[2] -= p2_t * epsilon_bar * delta_vec[1] * SVM_vec[1];

// Preschool booster (7th age group, age=5 yr, i = 6)
// need to change the equations for model version 4 still
if(add_preschool_booster) {
	DVM_vec[7] += p2_t * epsilon_bar * delta_vec[6] * MV_vec[6] + p2_t * epsilon_bar * delta_vec[6] * SVM_vec[6];
	DVMbar_vec[7] += p2_t * epsilon * delta_vec[6] * SVMbar_vec[6];
	DV_vec[7] += p2_t * epsilon * delta_vec[6] * SV_vec[6]; 
	DMV_vec[7] -= p2_t * epsilon_bar * delta_vec[6] * MV_vec[66];   
	DSV_vec[7] -= p2_t * epsilon * delta_vec[6] * SV_vec[6];
	DSVMbar_vec[7] -= p2_t * epsilon * delta_vec[6] * SVMbar_vec[6];
	DSVM_vec[7] -= p2_t * epsilon_bar * delta_vec[6] * SVM_vec[6];
}
//end_skel


//start_rsim
// Skeleton of stochastic model (=skeleton deterministic model)
int i, j, i2;

//Pointers to state variables
double *VM_vec = (double *) &VM_1;
double *VMbar_vec = (double *) &VMbar_1;
double *M_vec = (double *) &M_1;
double *MV_vec = (double *) &MV_1;
double *SMbar_vec = (double *) &SMbar_1;
double *SVMbar_vec = (double *) &SVMbar_1;
double *SVM_vec = (double *) &SVM_1;
double *S_vec = (double *) &S_1;
double *V_vec = (double *) &V_1;
double *SV_vec = (double *) &SV_1;
double *E_vec = (double *) &E_1;
double *I_vec = (double *) &I_1;
double *R_vec = (double *) &R_1;
double *CI1_vec = (double *) &CI1_1;
double *CI2_vec = (double *) &CI2_1;
double *CI3_vec = (double *) &CI3_1;
double *CI4_vec = (double *) &CI4_1;
double *CI5_vec = (double *) &CI5_1;

// Vector of parameters
double *q_vec = (double *) &q_1; // age-specific clinical fraction
double *N_vec = (double *) &N_1; // age-specific population size
double *delta_vec = (double *) &delta_1; // age-specific aging rates 
double lambda_vec[nages_mod];

// Calculate vaccine coverages
double p0_t = (t < tM) ? 0.0 : v0; // Maternal immunization
double p1_t = (t < tV) ? 0.0 : v1; // Primary series
double p2_t = (t < tV) ? 0.0 : v2; // Booster doses

// Re-create contact matrix
double CM_mat[nages_cmat][nages_cmat]; 
double *CM_vec = (double *) &CM_1; // Contact matrix stacked as vector (line by line), size nages_cmat*nages_cmat
for(i = 0; i < nages_cmat; i++) {
	for(j = 0; j < nages_cmat; j++) {
		CM_mat[i][j] = CM_vec[nages_cmat * i + j];
	}
}

// Compute seasonality factors in 5-10y and 10-20 y
CM_mat[1][1] *= seasC; // 5--10 y
CM_mat[2][2] *= seasT; // 10--15 y
CM_mat[3][3] *= seasT; // 15--20 y

// Create augmented contact matrix
double CM_mat_aug[nages_mod][nages_cmat]; // "Augmented" contact matrix with 2 extra rows for infant classes
double I_tilde[nages_cmat], N_tilde[nages_cmat]; // No infected and population size without infant age classes
I_tilde[0] = I_vec[0] + I_vec[1] + I_vec[2] + I_vec[3] + I_vec[4] + I_vec[5]; // Sum no infected in 0-4 mo, 4-12 mo, and 1-4 y
N_tilde[0] = N_vec[0] + N_vec[1] + N_vec[2] + N_vec[3] + N_vec[4] + N_vec[5]; // Sum population in 0-4 mo, 4-12 mo, and 1-4 y
for(i = 1; i < nages_cmat; i++) {
	I_tilde[i] = N_tilde[i] = 0.0;
	for(j = 1; j <= 5; j++) {
		I_tilde[i] += I_vec[5 * i + j];
		N_tilde[i] += N_vec[5 * i + j];
	}
}

// Add rows for the augmented contact matrix
for(i = 0; i < nages_mod; i++) {
	for(j = 0; j < nages_cmat; j++){
		if(i <= 5) {
			CM_mat_aug[i][j] = CM_mat[0][j];
		} else {
			i2 = (int) floor((i - 1) / 5); // Range: 1--14
			CM_mat_aug[i][j] = CM_mat[i2][j];
		}
	}
}

// Calculate force of infection in every age group
for(i = 0; i < nages_mod; i++) {
	lambda_vec[i] = 0.0; // Initialize values
	for(j = 0; j < nages_cmat; j++){
		lambda_vec[i] += CM_mat_aug[i][j] * (I_tilde[j] + iota) / N_tilde[j];
	}
	lambda_vec[i] *= 1 * q_vec[i]; // Adjust the R0 here with 0.5 and 1.5 for sensitivity analyses
}

// Blunting 
double epsilon_bar = (1.0 - b1) * epsilon; // Blunted vaccine effectiveness
double alphaV_bar = alphaV / (1.0 - b2); // Blunted duration of vaccine protection

if(debug) {
	Rprintf("t=%.2f, seasC = %.2f, seasT = %.2f\n\n", t, seasC, seasT);
}

// Transitions between compartments (one per compartment with the number of arrows that go out, start counting from 0)
// double from0_B[nages_mod][3];
double from0_M[nages_mod][4];
double from0_SMbar[nages_mod][4];
double from0_S[nages_mod][4];
double from0_VM[nages_mod][2];
double from0_MV[nages_mod][3];
double from0_VMbar[nages_mod][2];
double from0_SVMbar[nages_mod][3];
double from0_V[nages_mod][2];
double from0_SV[nages_mod][3];
double from0_SVM[nages_mod][3];
double from0_E[nages_mod][2];
double from0_I[nages_mod][2];
double from0_R[nages_mod][1];

// rates at age a
double rate1[1], trans1[1]; 
double rate2[2], trans2[2]; 
double rate3[3], trans3[3]; 
double rate4[4], trans4[4];


// Get three random birth rates 
double base_birth_rate = mu * N_tot * dt;
double BM_t = rpois(p0_t * epsilonM * base_birth_rate);
double BSMbar_t = rpois(p0_t * (1 - epsilonM) * base_birth_rate);
double BS_t = rpois((1 - p0_t) * base_birth_rate);

// Transition rates between compartments in newborns (age i = 0)
// B: second method to generate births, for now always gives consistent results with the three rpois approach
// rate3[0] = p0_t*epsilonM*mu;
// rate3[1] = p0_t*(1-epsilonM)*mu;
// rate3[2] = mu*(1-p0_t);
// reulermultinom(3, nearbyint(N_tot), &rate3[0], dt, &trans3[0]);
// for(j = 0; j<3; j++) from0_B[0][j] = trans3[j];
// M
rate4[0] = delta_vec[0]*(1-p1_t);
rate4[1] = tau;
rate4[2] = delta_vec[0]*p1_t*epsilon_bar;
rate4[3] = delta_vec[0]*p1_t*(1-epsilon_bar);
reulermultinom(4, nearbyint(M_vec[0]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_M[0][j] = trans4[j];
// SMbar
rate4[0] = delta_vec[0]*(1-p1_t);
rate4[1] = delta_vec[0]*p1_t*epsilon;
rate4[2] = delta_vec[0]*p1_t*(1-epsilon);
rate4[3] = lambda_vec[0];
reulermultinom(4, nearbyint(SMbar_vec[0]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_SMbar[0][j] = trans4[j];
// VM
rate2[0] = 0.0;
rate2[1] = 0.0;
reulermultinom(2, nearbyint(VM_vec[0]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_VM[0][j] = trans2[j];
// MV
rate3[0] = 0.0;
rate3[1] = 0.0;
rate3[2] = 0.0;
reulermultinom(3, nearbyint(MV_vec[0]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_MV[0][j] = trans3[j];
//VMbar
rate2[0] = 0.0;
rate2[1] = 0.0;
reulermultinom(2, nearbyint(VMbar_vec[0]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_VMbar[0][j] = trans2[j];
//SVMbar
rate3[0] = 0.0;
rate3[1] = 0.0;
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SVMbar_vec[0]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SVMbar[0][j] = trans3[j];
//SVM
rate3[0] = 0.0;
rate3[1] = 0.0;
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SVM_vec[0]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SVM[0][j] = trans3[j];
// S
rate4[0] = delta_vec[0]*(1-p1_t); //to S[1]
rate4[1] = delta_vec[0]*p1_t*epsilon; //to V[1]
rate4[2] = delta_vec[0]*p1_t*(1-epsilon); //to SV[1]
rate4[3] = lambda_vec[0];
reulermultinom(4, nearbyint(S_vec[0]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_S[0][j] = trans4[j];
//V
rate2[0] = 0.0;
rate2[1] = 0.0;
reulermultinom(2, nearbyint(V_vec[0]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_V[0][j] = trans2[j];
//SV
rate3[0] = 0.0;
rate3[1] = 0.0;
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SV_vec[0]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SV[0][j] = trans3[j];
//E
rate2[0] = delta_vec[0];
rate2[1] = sigma;
reulermultinom(2, nearbyint(E_vec[0]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_E[0][j] = trans2[j];
//I
rate2[0] = delta_vec[0];
rate2[1] = gamma;
reulermultinom(2, nearbyint(I_vec[0]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_I[0][j] = trans2[j];
//R
from0_R[0][0]=rbinom(nearbyint(R_vec[0]), pTrans(delta_vec[0], dt));
//if(t==0) {Rprintf("%.2f ", from0_R[0]);}



// Transition rates between compartments in the second age group (age i = 1, i.e., primary immunization)
// M
rate4[0] = delta_vec[1];
rate4[1] = tau;
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(M_vec[1]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_M[1][j] = trans4[j]; //not sure about this line
// SMbar
rate4[0] = delta_vec[1];
rate4[1] = lambda_vec[1];
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(SMbar_vec[1]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_SMbar[1][j] = trans4[j];
// VM
rate2[0] = delta_vec[1];
rate2[1] = alphaV_bar;
reulermultinom(2, nearbyint(VM_vec[1]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_VM[1][j] = trans2[j];
// MV
rate3[0] = delta_vec[1]*(1-p2_t*epsilon_bar);
rate3[1] = tau;
rate3[2] = delta_vec[1]*p2_t*epsilon_bar;
reulermultinom(3, nearbyint(MV_vec[1]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_MV[1][j] = trans3[j];
//VMbar
rate2[0] = delta_vec[1];
rate2[1] = alphaV;
reulermultinom(2, nearbyint(VMbar_vec[1]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_VMbar[1][j] = trans2[j];
//SVMbar
rate3[0] = delta_vec[1]*(1-p2_t*epsilon);
rate3[1] = delta_vec[1]*p2_t*epsilon;
rate3[2] = lambda_vec[1];
reulermultinom(3, nearbyint(SVMbar_vec[1]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SVMbar[1][j] = trans3[j];
//SVM
rate3[0] = delta_vec[1]*(1-p2_t*epsilon_bar);
rate3[1] = delta_vec[1]*p2_t*epsilon_bar;
rate3[2] = lambda_vec[1];
reulermultinom(3, nearbyint(SVM_vec[1]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SVM[1][j] = trans3[j];
// S
rate4[0] = delta_vec[1];
rate4[1] = lambda_vec[1];
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(S_vec[1]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_S[1][j] = trans4[j];
//V
rate2[0] = delta_vec[1];
rate2[1] = alphaV;
reulermultinom(2, nearbyint(V_vec[1]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_V[1][j] = trans2[j];
//SV
rate3[0] = delta_vec[1]*(1-p2_t*epsilon);
rate3[1] = delta_vec[1]*p2_t*epsilon;
rate3[2] = lambda_vec[1];
reulermultinom(3, nearbyint(SV_vec[1]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SV[1][j] = trans3[j];
//E
rate2[0] = delta_vec[1];
rate2[1] = sigma;
reulermultinom(2, nearbyint(E_vec[1]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_E[1][j] = trans2[j];
//I
rate2[0] = delta_vec[1];
rate2[1] = gamma;
reulermultinom(2, nearbyint(I_vec[1]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_I[1][j] = trans2[j];
//R
from0_R[1][0]=rbinom(nearbyint(R_vec[1]), pTrans(delta_vec[1], dt));


// Equations in the third age group (age i = 2, i.e., booster)
// Transition rates between compartments
// M
rate4[0] = delta_vec[2];
rate4[1] = tau;
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(M_vec[2]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_M[2][j] = trans4[j]; //not sure about this line
// SMbar
rate4[0] = delta_vec[2];
rate4[1] = lambda_vec[2];
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(SMbar_vec[2]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_SMbar[2][j] = trans4[j];
// VM
rate2[0] = delta_vec[2];
rate2[1] = alphaV_bar;
reulermultinom(2, nearbyint(VM_vec[2]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_VM[2][j] = trans2[j];
// MV
rate3[0] = delta_vec[2];
rate3[1] = tau;
rate3[2] = 0.0;
reulermultinom(3, nearbyint(MV_vec[2]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_MV[2][j] = trans3[j];
//VMbar
rate2[0] = delta_vec[2];
rate2[1] = alphaV;
reulermultinom(2, nearbyint(VMbar_vec[2]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_VMbar[2][j] = trans2[j];
//SVMbar
rate3[0] = delta_vec[2];
rate3[1] = lambda_vec[2];
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SVMbar_vec[2]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SVMbar[2][j] = trans3[j];
//SVM
rate3[0] = delta_vec[2];
rate3[1] = lambda_vec[2];
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SVM_vec[2]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SVM[2][j] = trans3[j];
//S
rate4[0] = delta_vec[2];
rate4[1] = lambda_vec[2];
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(S_vec[2]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_S[2][j] = trans4[j];
//V
rate2[0] = delta_vec[2];
rate2[1] = alphaV;
reulermultinom(2, nearbyint(V_vec[2]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_V[2][j] = trans2[j];
//SV
rate3[0] = delta_vec[2];
rate3[1] = lambda_vec[2];
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SV_vec[2]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SV[2][j] = trans3[j];
//E
rate2[0] = delta_vec[2];
rate2[1] = sigma;
reulermultinom(2, nearbyint(E_vec[2]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_E[2][j] = trans2[j];
//I
rate2[0] = delta_vec[2];
rate2[1] = gamma;
reulermultinom(2, nearbyint(I_vec[2]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_I[2][j] = trans2[j];
//R
from0_R[2][0]=rbinom(nearbyint(R_vec[2]), pTrans(delta_vec[2], dt));
//if(t==0) {Rprintf("%.2f ", from0_R[0]);}


// For all older age groups
for(i = 3; i < nages_mod; i++){
// Transitions rates between compartments
// M
rate4[0] = delta_vec[i];
rate4[1] = tau;
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(M_vec[i]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_M[i][j] = trans4[j]; //not sure about this line
// SMbar
rate4[0] = delta_vec[i];
rate4[1] = lambda_vec[i];
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(SMbar_vec[i]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_SMbar[i][j] = trans4[j];
// VM
rate2[0] = delta_vec[i];
rate2[1] = alphaV_bar;
reulermultinom(2, nearbyint(VM_vec[i]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_VM[i][j] = trans2[j];
// MV
rate3[0] = delta_vec[i];
rate3[1] = tau;
rate3[2] = 0.0;
reulermultinom(3, nearbyint(MV_vec[i]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_MV[i][j] = trans3[j];
//VMbar
rate2[0] = delta_vec[i];
rate2[1] = alphaV;
reulermultinom(2, nearbyint(VMbar_vec[i]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_VMbar[i][j] = trans2[j];
//SVMbar
rate3[0] = delta_vec[i];
rate3[1] = lambda_vec[i];
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SVMbar_vec[i]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SVMbar[i][j] = trans3[j];
//SVM
rate3[0] = delta_vec[i];
rate3[1] = lambda_vec[i];
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SVM_vec[i]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SVM[i][j] = trans3[j];
//S
rate4[0] = delta_vec[i];
rate4[1] = lambda_vec[i];
rate4[2] = 0.0;
rate4[3] = 0.0;
reulermultinom(4, nearbyint(S_vec[i]), &rate4[0], dt, &trans4[0]);
for(j = 0; j<4; j++) from0_S[i][j] = trans4[j];
//V
rate2[0] = delta_vec[i];
rate2[1] = alphaV;
reulermultinom(2, nearbyint(V_vec[i]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_V[i][j] = trans2[j];
//SV
rate3[0] = delta_vec[i];
rate3[1] = lambda_vec[i];
rate3[2] = 0.0;
reulermultinom(3, nearbyint(SV_vec[i]), &rate3[0], dt, &trans3[0]);
for(j = 0; j<3; j++) from0_SV[i][j] = trans3[j];
//E
rate2[0] = delta_vec[i];
rate2[1] = sigma;
reulermultinom(2, nearbyint(E_vec[i]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_E[i][j] = trans2[j];
//I
rate2[0] = delta_vec[i];
rate2[1] = gamma;
reulermultinom(2, nearbyint(I_vec[i]), &rate2[0], dt, &trans2[0]);
for(j = 0; j<2; j++) from0_I[i][j] = trans2[j];
//R
from0_R[i][0]=rbinom(nearbyint(R_vec[i]), pTrans(delta_vec[i], dt));
}

//Balance equations for the newborns (age i = 0)
VM_vec[0] += 0.0;
VMbar_vec[0] += 0.0;
M_vec[0] += BM_t - from0_M[0][0] - from0_M[0][1] - from0_M[0][2] - from0_M[0][3]; //from0_B[0][0] or BM_t
MV_vec[0] += 0.0;
SMbar_vec[0] += BSMbar_t + from0_M[0][1] - from0_SMbar[0][0] - from0_SMbar[0][1] - from0_SMbar[0][2] - from0_SMbar[0][3]; //from0_B[0][1] or BSMbar_t
SVMbar_vec[0] += 0.0;
SVM_vec[0] += 0.0;
S_vec[0] += BS_t - from0_S[0][0] - from0_S[0][1] - from0_S[0][2] - from0_S[0][3]; //from0_B[0][2] or BS_t 
V_vec[0] += 0.0;
SV_vec[0] += 0.0;
E_vec[0] += from0_SMbar[0][3] + from0_S[0][3] - from0_E[0][0] - from0_E[0][1]; 
I_vec[0] += from0_E[0][1] -  from0_I[0][0] - from0_I[0][1];
R_vec[0] += from0_I[0][1] - from0_R[0][0];
CI1_vec[0] += from0_SMbar[0][3];
CI2_vec[0] += from0_S[0][3];
CI3_vec[0] += 0.0;
CI4_vec[0] += 0.0;
CI5_vec[0] += 0.0;

//Balance equations for the second age group (age 1, i.e. primary immunization)
VM_vec[1] += from0_M[0][2] - from0_VM[1][0] - from0_VM[1][1]; 
VMbar_vec[1] += from0_SMbar[0][1] - from0_VMbar[1][0] - from0_VMbar[1][1];
M_vec[1] += from0_M[0][0] - from0_M[1][0] - from0_M[1][1];
MV_vec[1] += from0_M[0][3] - from0_MV[1][0] - from0_MV[1][1] - from0_MV[1][2];
SMbar_vec[1] += from0_SMbar[0][0] + from0_M[1][1] - from0_SMbar[1][0] - from0_SMbar[1][1];
SVMbar_vec[1] += from0_SMbar[0][2] + from0_VMbar[1][1] - from0_SVMbar[1][0] - from0_SVMbar[1][1] - from0_SVMbar[1][2];
SVM_vec[1] += from0_VM[1][1] + from0_MV[1][1] - from0_SVM[1][0] - from0_SVM[1][1] - from0_SVM[1][2]; 
S_vec[1] += from0_S[0][0] - from0_S[1][0] - from0_S[1][1];
V_vec[1] += from0_S[0][1] - from0_V[1][0] - from0_V[1][1];
SV_vec[1] += from0_S[0][2] + from0_V[1][1]- from0_SV[1][0] - from0_SV[1][1] - from0_SV[1][2];
E_vec[1] += from0_E[0][0] + from0_SVM[1][2] + from0_SV[1][2] + from0_SVMbar[1][2] + from0_S[1][1] + from0_SMbar[1][1]  - from0_E[1][0] - from0_E[1][1]; 
I_vec[1] += from0_I[0][0] + from0_E[1][1] - from0_I[1][0] - from0_I[1][1]; 
R_vec[1] += from0_R[0][0] + from0_I[1][1] - from0_R[1][0];
CI1_vec[1] += from0_SMbar[1][1];
CI2_vec[1] += from0_S[1][1];
CI3_vec[1] += from0_SV[1][2]; 
CI4_vec[1] += from0_SVMbar[1][2];  
CI5_vec[1] += from0_SVM[1][2]; 

//Balance equations for the third age group
VM_vec[2] += from0_SVM[1][1] + from0_MV[1][2] + from0_VM[1][0]- from0_VM[2][0] - from0_VM[2][1]; 
VMbar_vec[2] += from0_SVMbar[1][1] + from0_VMbar[1][0] - from0_VMbar[2][0] - from0_VMbar[2][1]; 
M_vec[2] += from0_M[1][0] - from0_M[2][0] - from0_M[2][1]; 
MV_vec[2] += from0_MV[1][0] - from0_MV[2][0] - from0_MV[2][1]; 
SMbar_vec[2] += from0_SMbar[1][0] + from0_M[2][1] - from0_SMbar[2][0] - from0_SMbar[2][1]; 
SVMbar_vec[2] += from0_SVMbar[1][0] + from0_VMbar[2][1] - from0_SVMbar[2][0] - from0_SVMbar[2][1]; 
SVM_vec[2] += from0_SVM[1][0] + from0_VM[2][1] + from0_MV[2][1] - from0_SVM[2][0] - from0_SVM[2][1]; 
S_vec[2] += from0_S[1][0] - from0_S[2][0] - from0_S[2][1]; 
V_vec[2] += from0_SV[1][1] + from0_V[1][0] - from0_V[2][0] - from0_V[2][1]; 
SV_vec[2] += from0_SV[1][0] + from0_V[2][1] - from0_SV[2][0] - from0_SV[2][1]; 
E_vec[2] += from0_E[1][0] + from0_SVM[2][1] + from0_SV[2][1] + from0_SVMbar[2][1] + from0_S[2][1] + from0_SMbar[2][1] - from0_E[2][0] - from0_E[2][1]; 
I_vec[2] += from0_I[1][0] + from0_E[2][1] -  from0_I[2][0] -  from0_I[2][1]; 
R_vec[2] += from0_R[1][0] + from0_I[2][1] - from0_R[2][0]; 
CI1_vec[2] += from0_SMbar[2][1]; 
CI2_vec[2] += from0_S[2][1]; 
CI3_vec[2] += from0_SV[2][1]; 
CI4_vec[2] += from0_SVMbar[2][1];  
CI5_vec[2] += from0_SVM[2][1]; 

for(i = 3; i < nages_mod; i++){
//Balance equations for the older age groups (i>4)
VM_vec[i] += from0_VM[i-1][0] - from0_VM[i][0] - from0_VM[i][1]; 
VMbar_vec[i] += from0_VMbar[i-1][0] - from0_VMbar[i][0] - from0_VMbar[i][1]; 
M_vec[i] += from0_M[i-1][0] - from0_M[i][0] - from0_M[i][1]; 
MV_vec[i] += from0_MV[i-1][0] - from0_MV[i][0] - from0_MV[i][1]; 
SMbar_vec[i] += from0_SMbar[i-1][0] + from0_M[i][1] - from0_SMbar[i][0] - from0_SMbar[i][1]; 
SVMbar_vec[i] += from0_SVMbar[i-1][0] + from0_VMbar[i][1] - from0_SVMbar[i][0] - from0_SVMbar[i][1]; 
SVM_vec[i] += from0_SVM[i-1][0] + from0_VM[i][1] + from0_MV[i][1] - from0_SVM[i][0] - from0_SVM[i][1]; 
S_vec[i] += from0_S[i-1][0] - from0_S[i][0] - from0_S[i][1]; 
V_vec[i] += from0_V[i-1][0] - from0_V[i][0] - from0_V[i][1]; 
SV_vec[i] += from0_SV[i-1][0] + from0_V[i][1] - from0_SV[i][0] - from0_SV[i][1]; 
E_vec[i] += from0_SVM[i][1] + from0_SV[i][1] + from0_SVMbar[i][1] + from0_S[i][1] + from0_SMbar[i][1] + from0_E[i-1][0] - from0_E[i][0] - from0_E[i][1]; 
I_vec[i] += from0_E[i][1] + from0_I[i-1][0] -  from0_I[i][0] - from0_I[i][1]; 
R_vec[i] +=  from0_R[i-1][0] + from0_I[i][1] - from0_R[i][0]; 
CI1_vec[i] += from0_SMbar[i][1]; 
CI2_vec[i] += from0_S[i][1]; 
CI3_vec[i] += from0_SV[i][1]; 
CI4_vec[i] += from0_SVMbar[i][1];  
CI5_vec[i] += from0_SVM[i][1]; 
}
//end_rsim